% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Edited by David Stern 17 March 2021 for final submission to JAERE.
% This code produces the Baseline Simulation in "Directed Technical Change
% and the British Industrial Revolution". 
% This version uses a normalized CES function and other equations are
% changed as required.
% This baseline version endogenously finds the two normalization constants.
% The counterfactual simulation uses the numerical values from this program
% as its normalization constants.

clc
close all
clear all 
tic

% Baseline Scenario
alpha=0.25; %share of energy in the production of the sectoral goods
beta=0.225; %sare of capital in the production of the sectoral goods
nu=(1-alpha-beta)/(1-beta); %within period decreasing returns to expenditure in R&D
p0=.4; %initial value of price ratio
Ns0=1; %initial number of varieties of Solow machines
Embar=1; %use of wood, constant over time
L0=1; %Initial population
T=18; %number of periods
theta=0; % Setting instead theta=beta^(beta/(beta-1))-1 allows simulation of 
% a model with competitive production of old varieties after the initial patents
% have expired. We used this version in earlier drafts of
% the paper. For the model presented in the paper Bm=Nm and Bs=Ns. More
% generally: Bm(t)=Nmt(t+1)+theta*Nmt(t) and Bs(t)=Nst(t+1)+theta*Nst(t).
sigma=4.0;%elasticity of substitution between sectoral goods
gamma=0.853;%share prameter in final good production
etas=2.26;%innovation productivity parameter in the Solow sector
etam=etas; %innovation productivity parameter in the Malthus sector
Nm0=5; %initialnumber of varieties of Malthus machines
esbar=0.115; %price of coal, constant over time
p0=gamma/(1-gamma); % the price ratio in the first period is predetermined in the normalized version
Embar=1; %use of wood, constant over time
Ns0=1; %initial number of varieties of Solow machines

%Import variable results from and export parameter values to function file
[yratio, Y,Ym,Ys,Ypc,pratio, Nmt, Nst, Bm, Bs, B, n, Em, Es, es, em, eratio, EI, SshrY, SshrL, SshrE, SshrEexp, gYpc,gNm,gNs,Lratio,L,ps,pm,Ls,Lm,Ymb,Ysb] = DTC_CD_Normalized_Baseline_fun(gamma,Nm0,etam,etas,sigma,beta,alpha,Ns0,Embar,esbar,nu,T,p0,theta,L0);

%Calibration algorithm. These equations compute the squared deviations from
%the stylized facts as discussed in Section 5.2. The model is fitted by
%varying the parameters Nm0, eta, esbar, and gamma and attempting to
%minimize crit.
e1=(log(em(1,1)/esbar)-log(2.))^2; 
e2=(log(Es(1,1)/Em(1,1))-log(0.3))^2; 
e3=(log(max(em)/em(1,1))-log(2.))^2; 
e4=(log(EI(T,1)/min(EI))-log(2.))^2;
e5=(log(Ypc(13,1)/Ypc(1,1))-log(2))^2;
e6=(log(Ypc(18,1)/Ypc(1,1))-log(5.4))^2;
crit=(e1+e2+e3+e4+e5+e6);

%Time variables are in 20 year steps, with main index, Year, starting in
%1540+20=1560. Year1 starts in 1520+20=1540 and is for Nm and Ns, for which
%t+1 means t for all other variables except growYpc. Year2 is for growYpc.
Year=1540+(1:T).*20;
Year1=1520+(1:(T+1)).*20;
Year2=1560+(1:(T-1)).*20;
%Half is used just to plot horizontal 0.5 line in Fig 9
Half=0.5*ones(1,T);
 
   
figure(1)
plot(Year,Ypc/Ypc(1),'linewidth',2);
title('CD model fig. 1')
xlabel('year')
ylabel('GDP per capita (1560=1)','FontSize',12)
axis([1550 1900 0 7])
 
figure(2)
plot(Year,em,'r:',Year,es,'r--',Year,eratio,'r-.','linewidth',2);
title('CD model fig. 2')
xlabel('year');
ylabel('energy price','FontSize',12);
legend('wood price','coal price','wood/coal price','location','northwest')
axis([1550 1900 0 5])
 
figure(3)
plot(Year,Em,'r:',Year,Es,'r--','linewidth',2);
title('CD model fig. 3')
xlabel('year');
ylabel('energy quantity','FontSize',12);
legend('wood quantity','coal quantity','location','northwest')
axis([1550 1900 0 70])
 
figure(4) % not crucial for comparing CD & LS?
plot(Year2,gYpc,'b-.','linewidth',2);
title('CD model fig. 4')
xlabel('year');
ylabel('growth rate of GDP per capita (in %/yr)','FontSize',12);
axis([1550 1900 -0.1 3])
 
figure(5)
plot(Year1,Nmt,'c:',Year1,Nst,'c-.','linewidth',2);
title('CD model fig. 5')
xlabel('year')
ylabel('number of varieties','FontSize',12)
legend('Malthus sector','Solow sector','location','northwest')
axis([1550 1900 0 30])
 
figure(7)
plot(Year,EI/EI(1),'r-.','linewidth',2);
title('CD model fig. 7')
xlabel('year');
ylabel('energy intensity (1560=1)','FontSize',12);  
axis([1550 1900 0 2])
 
figure(8)
plot(Year,yratio,'r-.','linewidth',2);
title('CD model fig. 8')
xlabel('year');
ylabel('Malthus/Solow output ratio, y','FontSize',12);  
axis([1550 1900 0 30])
 
figure(9) % was figure (17) in LS_JP15_30Dec14
plot(Year,SshrE,'r--',Year,SshrL,'g--',Year,SshrY,'k--',Year,SshrEexp,'r-.',Year,Half,'k:','linewidth',2);
title('CD model fig. 9')
xlabel('year');
ylabel('Solow sector shares','FontSize',12);
legend('S-share of energy','S-share of labor','S-share of output','S-share of energy cost','location','northwest');
axis([1550 1900 0 1])
 


% No figure (10) for CD as within-sector energy cost shares are constant